%  Copyright (C) 2002 Regents of the University of Michigan, portions used with permission 
%  For more information, see http://csem.engin.umich.edu/tools/swmf
\chapter{Organization and Modification of the Source Code \label{chapter:source_code}}

\section{\BATSRUS source \label{section:models}}

The source code of \BATSRUS\ is located in the {\tt src/},
Within the directory, the source code is split into
a number of different files that are named according to
their function and the the subroutines which they contain.

\section{Modules \label{section:modules}}

A module is a Fortran 90 programming unit that allows
definitions of global variables and data structures.  
While the module has other features, \BATSRUS\ mainly uses
modules to declare globally used variables.  For example the module 
{\tt ModMain} is contained in
the file {\tt ModMain.f90}.  We include parts of the module here
(...) indicates that we have removed lines) to show how modules are
used in \BATSRUS.
\begin{verbatim}
module ModMain

  character (len=10) :: CodeVersion='6.01'

  !\
  ! Block parameters.
  !/
  integer, parameter :: nBLK=150
  integer :: nMultiBlkLevels, globalBLK

  !\
  ! Cell array parameters.
  !/
  integer, parameter :: nCellsI=4,  &
                        nCellsJ=4,  &
                        nCellsK=4
  integer, parameter :: nI=nCellsI, nJ=nCellsJ, nK=nCellsK
  integer, parameter, dimension(3) :: nCells=(/nCellsI,nCellsJ,nCellsK/)
  integer, parameter :: gcn=2 ! number of ghostcells.
...
  !\
  ! Face array parameters.
  !/
  integer, parameter :: nFacesIx=nCellsI+1, &
                        nFacesJx=nCellsJ,   &
                        nFacesKx=nCellsK
  integer, parameter :: nFacesIy=nCellsI,   &
                        nFacesJy=nCellsJ+1, &
                        nFacesKy=nCellsK
  integer, parameter :: nFacesIz=nCellsI,   &
                        nFacesJz=nCellsJ,   &
                        nFacesKz=nCellsK+1
...
  !\
  ! Geometry parameters.
  !/
  real  ::    x1, x2, y1, y2, z1, z2
...
  !\
  ! Block cell-centered MHD solution array definitions.
  !/
  real,  dimension(1-gcn:nCellsI+gcn, &
                   1-gcn:nCellsJ+gcn, &
                   1-gcn:nCellsK+gcn,nBLK) :: &
    x_BLK,y_BLK,z_BLK,R_BLK,R2_BLK, &           !Block cell coordinates
    rho_BLK,rhoUx_BLK,rhoUy_BLK,rhoUz_BLK, &    !Block cell state
    Bx_BLK,By_BLK,Bz_BLK,E_BLK,p_BLK, &         
    B0xCell_BLK,B0yCell_BLK,B0zCell_BLK, &      !Block intrinsic field
    time_BLK, &                                 !Block time 
    tmp1_BLK, tmp2_BLK                          !Temporary block variables
...
  !\
  ! Other block solution and geometry parameters defined for each block.
  !   cell dimensions, minimum radius of a block, face areas, cell volume and 
  !   inverse of the cell volume.  
  !/
  real, dimension(nBLK) :: dx_BLK, dy_BLK, dz_BLK, Rmin_BLK, Rmin2_BLK, &
			    fAx_BLK, fAy_BLK, fAz_BLK, cV_BLK, cVinv_BLK

  integer, dimension(nBLK) :: global_block_number

  !\
  ! X Face local MHD solution array definitions.
  !/
  real, dimension(0:nFacesIx+1,0:nFacesJx+1,0:nFacesKx+1) ::     &
    rhoFaceL_x,UxFaceL_x,UyFaceL_x,UzFaceL_x,              & !\
    BxFaceL_x,ByFaceL_x,BzFaceL_x,pFaceL_x,                & ! Face Left X
    rhoFaceR_x,UxFaceR_x,UyFaceR_x,UzFaceR_x,              & !\
    BxFaceR_x,ByFaceR_x,BzFaceR_x,pFaceR_x,                & ! Face Right X
    rhoFaceF_x,rhoUxFaceF_x,rhoUyFaceF_x,rhoUzFaceF_x,     & !\
    BxFaceF_x,ByFaceF_x,BzFaceF_x,EFaceF_x,                & ! Face Flux X
    VdtFace_XI                                                !V/dt Face X
...
  ! Logical for bodies
  logical :: body1=.false., body2=.false.
...
end module ModMain
\end{verbatim}
While much of {\tt ModMain} has been left out, what remains above
shows the typical way in which modules are used in \BATSRUS.  {\tt ModMain}
defines the main variable structure of \BATSRUS\ that contains the
blocks in which the solution is stored. It additionally defines any
variable that has an effect on computation and therefore must
be accessable to the majority of the code.  

The variables {\tt nBLK} and {\tt nCells} define the number of solution
blocks and the size of blocks it terms of grid cells.  
These variables are then used to define blocks for each of
the solution variables ({\tt x\_BLK, rho\_BLK, rhoUx\_BLK, Bx\_BLK, p\_BLK} 
for example) and for face fluxes ({\tt rhoFaceL\_x, UxFaceL\_x, pFaceL\_x}).
In addition, variables like {\tt body1} and {\tt body2} are necessary so that
the code knows if there are bodies at which it should apply boundary
conditions.

In addition to {\tt ModMain}, there are other modules.  Each of these contains
variable definitions for related quantities.  For example, quantities 
associated
with input and output are found in {\tt ModIO}, ionospheric quantities are
found in {\tt ModIonosphere} and physical quantities like fundamental
constants are found in {\tt ModPhysics}.

Any subroutine that needs to know about the globally defined 
variables must include the module using the {\tt use}
command
\begin{verbatim}
use ModMain
\end{verbatim}

In addition to the above description, modules can be used to define
what is called an interface.  This allows a user to call a single 
routine which will behave slightly differently when called
with different sets of arguments.  There are a few instances of this
use in {\tt ModInterface}.

\section{Parallel Programming and Blocks \label{section:parallel_programming}}

In order to understand the source code of \BATSRUS\
it is important to understand how the code runs on
multiple processors and how the block structure influences
coding methods.  The first principle to understand is
that the code runs independently on each processor.  All memory
is local to the processor and the only way that information is
shared between processors is through the passing of messages.
What this means is that the same variables are declared on each
processor and are independently loaded.  For example, in 
{\tt ModMain} above, the variable {\tt nBLK} is the number of
blocks that are stored in memory on this processor.  On another
processor {\tt nBLK} is also defined and set to the same value.
Variables that are declared in modules are global variables 
accesable to any routine that ``uses'' them, however, being
global means global to a specific processor.  For example,
the variable {\tt globalBLK} is the block number of the 
block being processed currently on this processor.  In contrast,
the variable {\tt global\_block\_number} holds a unique number for
each block across all processors.  To
load this value, blocks have to recieve information from other
processors.  This aspect of parallel programming means that 
most of the time code can be writen as if it were intended
for use on a single processor.  Only when information needs to
be shared between processors (like
{\tt global\_block\_number} definitions or ghost cell values) does 
communication and parallel message passing come into play.
Many changes can be made without having to 
delve deeply into MPI and message passing.

The second important principle to understand is that \BATSRUS\
is built around blocks.  While this fact seems trivial
and has now been discussed several times (see section~\ref{section:grid}),
it is fundamental to the coding of \BATSRUS.  Because of the block 
based structure most routines are designed to perform an algorithm
on the simple Cartesian mesh of a single block.  Once an algorith is
 working on a single block, a loop over all the blocks is added and the
procedure now operates on the entire grid.  The structured
Cartesian mesh is easier to program for than an unstructured
mesh and the loop over blocks takes care of the rest.  The
fact that blocks are different size doesn't matter since the
general structure of each block is the same.

Finally, it is important to understand that since the program
is executing independently on each processor  output routines
and write statements are executed on each processor.  In the
code, the variable {\tt me\_world} contains the processor number
of each processor.  The test
\begin{verbatim}
if (me_world == 0) then 
...
end if
\end{verbatim}
is often used when writing to the screen or to a log file so
that only one entry (not one per processor) is made.

\section{Some Specific Global Variables \label{section:variables}}

In order to understand many of the routines in \BATSRUS\ it is helpful
to know  the meaning of some of the global variables.  In general,
variables that are named similarly are related and are dimensioned the
same.  There are, of course, exceptions to this rule.  In the following,
if a variable {\tt x\_BLK} is explained, the user can assume that
the varibles {\tt y\_BLK} and {\tt z\_BLK} are similarly defined.
Consult the module files to find out the names of the other similarly defined solution
variables.  Fortran 90 is not case sensitive, but many variables are writen
with capital letters to make them easy to read.  {\bf Be aware that there may be 
some inconsistancies in the use of capitalization in different parts of the code.}

\begin{tabbing}
{\bf numprocs} \hspace{1.0in} \= Number of processors being used in this run. \\
{\bf me\_world}               \> Unique number of each processor. Numbered from 0 to numprocs-1 \\
{\bf Srho, SrhoUx, ...}       \> Source terms \\
{\bf gradX\_rho, gradX\_Ux, ...} \> Gradients in the X, Y and Z directions \\
{\bf B0...}                   \> Arrays storing the intrinsic magnetic field of the central body \\
{\bf B...}                    \> Array storing B$_1$, the deviative part of the magnetic field (see
                                 section~\ref{section:B0}). \\
{\bf nBLK}                    \> Maximum number of blocks per processor. \\
{\bf nMultiBlockLevels}       \> The maximum number of used blocks on any of the processors. \\
{\bf nCellsI,nCellsJ, nCellsK}\> Number of computational cells in the i,j and k directions. \\
{\bf nI,nJ, nK}               \> Same as {\tt nCellsI}, {\tt nCellsJ} and {\tt nCellsk}. \\
{\bf nFacesIX, ..., nFacesIy} \> Number of x faces in the i direction, ..., \\
                              \> Number of y faces in the i direction, ... (see below). 
\end{tabbing}

\noindent
{\bf gcn} \\
The number of ghost cells on each face. This variable is set to 2.
This means that the actual dimensions of a block is 
{\tt (nCellsI+2*gcn)x(nCellsj+2*gcn)x(nCellsk+2*gcn)}.
The variable {\tt gcn} is used inconsistently throughout the code. 
Sometimes {\tt gcn} is used, but other times the value 2 is hard coded.
Because of this the value of gcn should not be changed.
 
\ \ \\
{\bf globalBLK, iBLK}.  \\
These variables are for looping over all the blocks on a 
single processor.  They contains the number of the block currently being processed 
on this processor.  The variable {\tt globalBLK} is global and therefore
does not have to be passed to routines that use {\tt ModMain}, the variable 
{\tt iBLK} is not global and is passed into routines.

\ \ \\
{\bf \_BLK, BLK} \\
Variables that end with {\tt \_BLK} or {\tt BLK} are dimensioned to hold
a value for each block on that processor.  Remember that {\tt nBLK} is the
number of block per processor.  The variables declaration above
in {\tt ModMain} include the following definition:
\begin{verbatim}
  real,  dimension(1-gcn:nCellsI+gcn, &
                   1-gcn:nCellsJ+gcn, &
                   1-gcn:nCellsK+gcn,nBLK) :: &
    x_BLK,y_BLK,z_BLK,R_BLK,R2_BLK, &           !Block cell coordinates
\end{verbatim}
Here the variable {\tt x\_BLK} is the x position of the center of each cell.
The first three dimensions are the dimensions of a single block, including
ghost cells.  The forth dimension indicates that this variable stores the
locations for every block on this processor.  Similarly,
\begin{verbatim}  
  real, dimension(nBLK) :: dx_BLK, dy_BLK, dz_BLK, Rmin_BLK, Rmin2_BLK, &
\end{verbatim}
defines the cell size, {\tt dx\_BLK}, for each block on this processor.  Since
the cell size is constant for each block this variable is dimensioned 
simply by the number of blocks allowed per processor, {\tt nBLK}.

\ \ \\
{\bf unusedBLK} \\
This is a logical variable that indicates whether or not a block is
currently being used or not.  It is dimensioned by {\tt nBLK} and is
{\tt .true.} when the block is unused.

\ \ \\
{\bf rhoFaceL\_x, rhoFaceR\_x, rhoFaceF\_x ...}  \\
The face variables, like 
{\tt rhoFaceL\_x} and {\tt rhoFaceR\_x} contain the values of rho at the
interface between two cells.  The ``L'' and ``R'' in the name indicate
the value at the face to the left (the direction of lower index) and
to the right (the direction of higher index) of the interface between two
cells. 
The ``F'' indicates that this is the flux through the face. The ``x''
indicates that these are faces whose normal vector is in the x direction.
The reason that there is a left and right value at each face 
is described in section~\ref{section:finite_volume}.  There can be
only one flux through any given face.  The variables are
are dimensioned according to
\begin{verbatim}
  integer, parameter :: nFacesIx=nCellsI+1, &
                        nFacesJx=nCellsJ,   &
                        nFacesKx=nCellsK
  integer, parameter :: nFacesIy=nCellsI,   &
                        nFacesJy=nCellsJ+1, &
                        nFacesKy=nCellsK
  integer, parameter :: nFacesIz=nCellsI,   &
                        nFacesJz=nCellsJ,   &
                        nFacesKz=nCellsK+1

  !X Faces
  real, dimension(0:nFacesIx+1,0:nFacesJx+1,0:nFacesKx+1) ::     &
    rhoFaceL_x, rhoFaceR_x, rhoFaceF_x

  !Y Faces
  real, dimension(0:nFacesIy+1,0:nFacesJy+1,0:nFacesKy+1) ::     &
    rhoFaceL_y, rhoFaceR_y, rhoFaceF_y

  !Z Faces
  real, dimension(0:nFacesIz+1,0:nFacesJz+1,0:nFacesKz+1) ::     &
    rhoFaceL_z, rhoFaceR_z, rhoFaceF_z
\end{verbatim}
Notice that there are a different number of x faces in the x direction
than in the y and z directions. Figures~\ref{fig:faces1} and \ref{fig:faces2}
show how cell centers and faces are numbered.
\begin{figure}
\begin{center}
\includegraphics*[width=2.5in]{faces.pdf}
\end{center}
\caption{2D numbering of faces.  A single block is shown with nCellsI=4, nCellsJ=4.
The computational cells are those inside the heavy black line.  The ghost cells are those
outside the line.  The X faces for which {\tt FaceL\_x} and {\tt FaceR\_x} are defined
are indicated by the thick gray line.  Y faces and Z faces
are defined similarly.  The naming convention for the left and right faces is given by 
L$_{i,j}$ and R$_{i,j}$ respectively.}
\label{fig:faces1}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics*[width=2.5in]{face_numbering.pdf}
\end{center}
\caption{The 2D numbering of faces relative to cell centers 
is shown for a single cell and its 
neghbors.  Notice that the labeling of ``L'' and ``R'' is relative to the face not
the cell center.}  
\label{fig:faces2}
\end{figure}

\ \ \\
{\bf Neighbors} \\
In \BATSRUS\ neigbors refer to blocks that are bordering the current
block.  It is important to know which block is a neighboring block
and what the refinement level of that block is so that messages and
ghost cell information can be passed correctly. All neighbor
variables are dimensioned by {\tt nBLK} and tell you the properties
of a neighboring cell in a certain direction.  The neighbor directions
in \BATSRUS\ are defined according to:
\begin{tabbing}
{\bf east}  \hspace{.25in}\= In the minimum X direction \\
{\bf west}                \> In the maximum X direction \\
{\bf south}               \> In the minimum Y direction \\
{\bf north}               \> In the maximum Y direction \\
{\bf bot}                 \> In the minimum Z direction \\
{\bf top}                 \> In the maximum Z direction
\end{tabbing}

The properties of a neighboring cell that must be known are the
refiment level, the processor on which the neighbor is located,
the block number on that processor of the neighbor, and if the
child number of that neighbor.  We discuss
the level of neighbor first.  Blocks will smaller cells are higher
in level, meaning that they have been refined more.  The level
values of neighbors are
\begin{tabbing}
 0 \hspace{.5in}   \= Neighbor at same refinement level (same size cells). \\
-1                 \> Neighbor more refined (smaller cells). \\
+1                 \> Neighbor coarser (larger cells). \\
NOBLK              \> Neighbor does not exist or is not used.
\end{tabbing}

Neighbors may be found on another processor.  It is important to
know which processor and which block on the processor is the neighbor.
The variables for neighbor blocks store the block number on the
processor where the block is located.

There are currently three sets of redundant variables for 
neighbors in \BATSRUS.  These are slightly different and therefore are
usefull in different situations.  The set of variables that
contains information about neighbors in all directions including
corners is:
\begin{verbatim}
  integer, dimension( -1:1, -1:1, -1:1, 4, nBLK) :: &
     BLKneighborPE, BLKneighborBLK, BLKneighborCHILD
  integer, dimension( -1:1, -1:1, -1:1, nBLK) :: BLKneighborLEV
\end{verbatim}
The dimensions for  {\tt BLKneighborPE}, {\tt BLKneighborBLK}
and  {\tt BLKneighborCHILD}
correspond to the following: (i,j,k,\#\_of\_neighbors,block). Here i,j,k are
indexes to look up or down in the x,y,z direction.  With this
structure it is possible to look in any of the 27 neighbor directions 
(faces, edges and corners)
just by choosing the indexes correctly.  
Since only one level change is permitted
between neighboring solution blocks, there are either 1 or 4 
neighboring blocks in each  direction. This is the
reason for the \#\_of\_neighbors dimension.
The {\tt  BLKneighborLEV} variable is not defined with this 
dimension since it tells only whether the block is refined
 or not without giving
specific information about location of the block.

The other two sets of variables differ from the previous in that
the contain information only about the 6 cardinal direction and 
no information about corners of edges.  The set differ from
each other only in compactness.  Meaning that the single array
form uses an index to contain all the individually defined variables.
The level variables are:
\begin{verbatim}
  integer, dimension(nBLK) :: &
       neiLtop, neiLbot, neiLeast, neiLwest, neiLnorth, neiLsouth
  integer, dimension(east_:top_,nBLK):: neiLEV
\end{verbatim}
The variables for processor (P, PE) and block (B, BLK) are:
\begin{verbatim}
  integer, dimension(4,nBLK) :: &
       neiPtop, neiPbot, neiPeast, neiPwest, neiPnorth, neiPsouth, &
       neiBtop, neiBbot, neiBeast, neiBwest, neiBnorth, neiBsouth
  integer, dimension(4,east_:top_,nBLK) :: neiPE, neiBLK
\end{verbatim}

\ \ \\
{\bf body\_BLK, true\_BLK, true\_cell} \\
These variables are used for testing whether or not a block or a cell
is special because it touches or is inside a body or touches the
outer boundary.  The variables are defined as follows:
\begin{tabbing}
body\_BLK \hspace{.25in} \= Defined for each block.  This variable is {\tt .true.} when any cell \\ 
                         \> \hspace{.25in}\= (including ghost cells) in the block is inside a body. \\
true\_cell               \> Defined for each cell in each block.  This variable is {\tt .true.} \\
                         \>     \> when a cell is not inside the body or beyond the outer
                                   boundary. \\    
true\_BLK                \> Defined for each Block.  This variable is {\tt .true.} when all cells \\
                         \>     \> (including ghost cells) in the block are true\_cells
\end{tabbing}

\section{Finding Default Settings \label{section:find_defaults}}

Many variables are set to default settings.  This allows the user
to run the code without having to set all the complicated values 
that must be set.  Eventually, these may have to be changed and 
then the problem of finding where they are set presents it self.

Many parameters are problem type dependent (for example, solar wind input
values).  These are set on a case by case basis in {\tt read\_inputs.f90}.
Other parameters are for controlling details of the numerics and generally
do not need to be changed.  Some of these are set in {\tt read\_inputs.f90}
while others are set in the various modules ({\tt ModMain} and {\tt ModProject} for
example).  Finally, physical constants and normalizations are set in
{\tt ModPhysics} and {\tt set\_physics.f90}.

\section{Normalization \label{section:normalization_source}}

Section~\ref{section:normalization} gave the detailed functions by
which \BATSRUS\ normalizes the solution variables.  In this section
we provide the user with detailed information about how normalization
is done and what changes the user might wish to make.  As mentioned in
section~\ref{section:normalization}, normalization of all quantities is
done using 
\begin{tabbing}
$\rho_0$  \hspace{0.5in} \= A reference density \\
$u_0$                    \> A reference velocity \\
$R_0$                    \> The typical length scale of the problem
\end{tabbing}
The values assigned to each of these quantities is problem type dependent.
Once these three variables have been assigned, the normalizations for all
other quantities can be computed.  The normalization factors and the initial
normalizations are done in routines {\tt set\_physics\_constants} and 
{\tt set\_dimensions} in the file {\tt set\_physics.f90}. 

{\bf There are two sets of variables for doing normalization}.  As an example we look
at {\tt unitSI\_B} and {unitUSER\_B} to represent the two sets.  One set of
variables ({\tt unitSI}) is for normalizing in standard SI units 
(m, m$^{-3}$, kg, s, m/s, T, V, Pa, J).  {\bf The {\tt unitSI} variables are consistant and
can be combine to get the normalization for any other SI variable}.  
The other set or variables ({\tt unitUSER}) contains normalization factors
in the users prefered units (for example: planetary radii, cm$^{-3}$, km/s, nT, mV, nPa).
{\bf This set of units may or may not be self consistant and care should be taken
when make combinations}.

The two sets of units are used to normalize input parameters and then to 
dimensionalize \BATSRUS\ variables for output.  Users generally prefer to
input parameters and create output in non-SI units.  With velocity, for example,
values of 400 km/s are typical in the solar wind.  Inputing and outputing
values in m/s is less convenient than in km/s.  The {\tt unitUSER} variables
are used for this purpose.  The desired quantities are load by problem type in
the routines {\tt set\_physics\_constants} and 
{\tt set\_dimensions} in the file {\tt set\_physics.f90}.  

\ \ \\
{\bf Normalizing input in the user's units for use in \BATSRUS\ takes the form:}
\begin{verbatim}
output of magnetic field in nT = unitUSER_B * B_BLK
\end{verbatim}
{\bf Dimensionalizing \BATSRUS\ variables for output in the user's units takes the form:}
\begin{verbatim}
B_BLK = input of magnetic field in nT / unitUSER_B
\end{verbatim}
where {\tt B\_BLK} is the variable in \BATSRUS\ that stores the unitless (normalized)
magnetic field.  The {\tt unitSI} variables are used similarly, but are in the
standard SI units.

\section{Testing, Debugging and {\tt set\_okTest} \label{section:tests}}

Controlling output while testing and debuggig can be challenging in \BATSRUS\ 
because of the multiple processors and multiple blocks.  
The testing input commands (see the USER MANUAL) can make getting useful 
information a little simpler.  The basis for controlling output is the 
command {\tt \#TEST} and the routine {\tt set\_oktest}.  As outlined in 
the USER MANUAL, this command takes the form
\begin{verbatim}
#TEST  
exchange_messages read_inputs
\end{verbatim}
The string ``{\tt exchange\_messages read\_inputs}'' is the {\tt test\_string}.
This line will be parsed using spaces so individual strings should be separated
by spaces.  The {\tt test\_string} can be combined with the routine 
{\tt set\_oktest} to control output.  A call to {\tt set\_oktest} looks like
\begin{verbatim}
  call set_oktest(``some_string'',oktest,oktest_me)
\end{verbatim}
where ``{\tt some\_string}'' is a string and {\tt oktest} and 
{\tt oktest\_me} are logicals which are set in {\tt set\_oktest} as follows:
\begin{verbatim}
  oktest=index(test_string,oktest_string)>0
  oktest_me = oktest .and. me_world==PROCtest
\end{verbatim}
This means that {\tt oktest} is {\tt .true.} if the 
{\tt oktest\_string} is contained in the {\tt test\_string} from the
{\tt \#TEST} command.  
The logical {\tt oktest\_me} is {\tt .true.} if 
{\tt oktest} is {\tt .true.} and {\tt me\_world == PROCtest}, or in other words,
if this processor has a processor number ({\tt me\_world}) equal to {\tt PROCtest}.

The usefulness of this routine is outlined in following example. 
If the routine {\tt exchange\_messages} contains the following
\begin{verbatim}
  logical :: oktest, oktest_me, oktime, oktime_me
  call set_oktest('exchange_messages',oktest,oktest_me)
  if(oktest)write(*,*)'Checked negative P, me=',me_world
\end{verbatim}
Then when the {\tt \#TEST} command is called with
\begin{verbatim}
#TEST  
exchange_messages            Test String
\end{verbatim}
\BATSRUS\ will the write the message ``Checked negatife P''.  If the {\tt test\_string}
does not contain the string ``exchange\_messages'' then no outpu will be writen.
The common usage for {\tt oktest} is to put the name of a subroutine 
in the {\tt test\_string} line.  This will generally cause output
to be written from inside that subroutine that is useful in debugging.
However, the {\tt test\_string} can contain any strings in any combination
that the user find usefull.


The routine {\tt set\_oktest} can be used in conjunction other testing commands to 
further customize output.  For example setting
\begin{verbatim}
#TEST  
correctP                 Test String

#TESTIJK
1                        Itest          (cell index for testing)
1                        Jtest          (cell index for testing)
1                        Ktest          (cell index for testing)
1                        BLKtest        (block index for testing)
0                        PROCtest       (processor index for testing)
\end{verbatim}
in the input file and then including the 
\begin{verbatim}
logical :: oktest, oktest_me

if(globalBLK==BLKtest)then
   call set_oktest('correctP',oktest,oktest_me)
else
   oktest=.false.; oktest_me=.false.
end if

f(oktest_me.and.i==Itest.and.J==Jtest.and.K==Ktest)&
   write(*,*)'CorrectP, initial P,E=',qp,qe
\end{verbatim}
in the {\tt correctP} routine will produce the ``CorrectP, initial P,E=''
error message only when {i,j,k=1,1,1} in block 1 on processor 0.  

\section{Location of Physics in the Source Code \label{section:locating_physics}}

The majority of files and routines that make up \BATSRUS\ are about numerics,
AMR and message passing.  The {\em physics} of \BATSRUS\, or in other words, the
routines and constants that define specific problems and model specific physical
processes can be found scattered over a few routines.  Without explaining
all the details, in this section we point out where certain physical parameters
are loaded or computed.

%\parskip=.25in \parindent=0pt

\ \ \\
{\bf ModPhysics.f90} \\
This module contains the definition of physical constants and variable declarations
for {\tt unitUSER} and {\tt unitSI} variables.  Other variable declarations are found
here for solar wind parameters, the intrinsic magnetic field, the problem specific
variables like {\tt R\_saturn}, the radius of Saturn.

\ \ \\
{\bf set\_physics.f90} \\
The routines contained in this file set physical constants that depend on
input parameters (like the speed of light - see section~\ref{section:boris})
and also sets the variables  {\tt unitUSER} and {\tt unitSI} for doing 
normalization of input and output.

\ \ \\
{\bf set\_ICs.f90 and create\_soln\_block.f90} \\
This routine sets the initial condition on the entire grid.  The initial setting
is problem dependent.

\ \ \\
{\bf set\_BCs.f90} \\
This routine sets the boundary conditions.  This routine is quite general and
most changes to boundary condition can be made in the input file.  However,
boundary conditions can drive the solution of a problem and changes often
need to be made to allow the user to apply a previously unused set of boundary
condtions.

\ \ \\
{\bf calc\_sources.f90} \\
This routine loads the source term arrays that are used to update the solution.
These arrays can be loaded with values that mimic physical processes like mass
loading (see section~\ref{section:mass_loading}). 

\ \ \\
{\bf calculate\_dipole\_tilt.f90} \\
This routine calculates the dipole and rotation axis of the central body in
terms of the date and time.  Specifics of the rotation and the axis of the intrinsic
magnetic field of the body are found here.

\ \ \\
{\bf B0multipole\_averages.f90} \\
This routine calculates cell averaged and face averaged values of the intrinsic
magnetic field.  While values to load the multipole expansion of the field can
be loaded in {\tt set\_physics.f90}, the user may need to change the way that
magnetic field are calculated for some problems.

\ \ \\ 
{\bf body\_force\_averages.f90} \\
This routine calculates cell averaged and face averaged values of the 
gravitational field of the central body and also calculates coriolis forces for
the rotating frame calculation of the heliosphere.

\ \ \\
{\bf specify\_refinement.f90 and amr\_criteria.f90}
These routines specify the initial refinement and the criteria that are
used to do automatic refinement (AMR).  Choosing initial refinement 
types and AMR criteria is
controlled from the {\tt PARAM.in} file and many options are available.
However, the user may which to change the initial grid structure or
add a new physical quantity to use for AMR.

\ \ \\
{\bf magnetosphere.f90 and srcIONO/} \\
These contain specifics about how the field aligned currents are mapped to the 
ionosphere and how the ionospheric potential is solved.
%\parskip=0pt \parindent=3em


\section{Adding Problem Types \label{section:problem_type}}

\BATSRUS\ was designed to be able to solve many different astrophysical problems.
The {\tt problem\_type} variable allows the code to contain the necessary 
modifications to perform simulations of many different problems.
The best way to add a new problem is to follow an existing problem to find out
what it does and in which routines specifics are contained.
Typing
\begin{verbatim}
grep problem_type *.f90 
\end{verbatim}
in the {\tt src/} directory will give a listing of all the occurances of the 
{\tt problem\_type} variable in the source code.  Looking at each of these
and adding the appropriate code for the new application is the best method
for solving a new problem with \BATSRUS.

When making changes, keeping the {\tt problem\_type} variable confined to 
the fewest routines possible will keep the code general and make adding problem
types earier.

\section{Adding Input Parameters \label{section:adding_input}}

Input parameters are read in the {\tt read\_inputs.f90} routine.  The format of
input entries is described in the USER MANUAL.  The best
way to add new input paramters is to look in {\tt read\_inputs.f90} and follow the many
examples.  The parameters should be declared in the appropriate module 
(e.~g. {\tt ModMain}) and default values should be set in the module or in 
{\tt set\_defaults} or {\tt set\_problem\_defaults} in {\tt read\_inputs.f90}

\section{Adding Source Files \label{section:adding_files}}

A new source file that needs to be compiled with \BATSRUS\ should be
put in the {\tt src/} directory.  The file will need to be added to the
list of file in the {\tt src/Makefile.COMMON} file.








